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ABSTRACT 


Numerical analysis of two dimensional steady 
state unconfined seepage through partially lined channels 
has been presented. The channels are either lined at base 
only and sides unlined or sides only lined and base unlined. 

The soil is considered to be homogenous and isotropic while 
the subsoil was either impermeable ( condition B ) or infi- 
nitely permeable ( condition A ). The parameters governing 
seepage are the depth of water in the channel » the base width 
of the channel, the position of water table and the position 
of subsoil layer. The side slope of the channel is kept cons- 
tant at 1H:1V. The condition A is a more active state of 
seepage than condition B. There is a marginal effect of 
position of water table on the percentage reduction in dis- 
charge for both cases of lining for condition B while for 
condition A the percentage reduction in seepage is found to 
increase with the lowering of the water table for base only 
lined case and decrease for sides only lined case. The 
lowering of the impermeable layer increased the percentage 
reduction in seepage for both cases of lining for condition B 
while lowering of the permeable layer decreased the percentage 
reduction in seepage for condition A. For less than 0.4 

for a set of other parameters the base lining only of the 
channel was found to be more effective in reducing seepage 
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while for greater than 0.4 sides lining was found to 

be more effective for condition A, It is also observed that 
the percentage reduction in seepage due to lining is not 
linearly related to the percentage of the total v/etted peri- 
rueter lined. The analysis is done with the help of a computer 
program based on finite element method. 
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CHAPTER 1 

INTRODUCTION 

1 . 1 General ; 

Canals are the main links between the storage and 
the user of the stored water. Since there is a large area 
under cultivation in our country importance of irrigation 
cannot be overemphasized. And indeed channel design and 
maintenance have gained considerable attention of investi- 
gators in this field. Though canals are essential for 
irrigation they also cause certain environmental disturbance 
by the waterlogging due to excessive seepage through their 
perimeter, 

A considerable quantity of water flowing in irriga- 
tion channels is lost to the adjacent areas through seepage. 

It has been estimated that transit losses in the alluvial 
channels in the northern India are 17 percent for the main 
canals and branches, 8 percent for distributaries, and 20 
percent for the watercourses, which gives a total loss of 
45 percent of the water entering the canal head. The nece- 
ssity of lining of irrigation channels cannot, therefore, be 
overemphasized. Besides reducing seepage losses in channels 
and thereby leading to more efficient utilization, the lining 
prevents waterlogging and weed growth, and reduces require- 
ment of land and maintenance cost. It has been estimated that 
lining of all the existing unlined channels in India could 



save enough v/ater to irrigate an additional six million 
hectares against the present irrigated area of 36,4 million 
hectares [38], 

However, nowaday s as a matter of policy all major 
channels are fully lined to eliminate the problem of v/ater- 
logging due to seepage and to improve the efficiency of the 
channels. Therefore the main cause of concern are the remain- 
ing unlined alluvial channels built in the beginning of this 
century which cause heavy seepage and accompanying water logg- 
ing. To mitigate the seventy of the problem it is necessary 
to partially line their perimeter, at least as a first phase, 
Partial lining may also be resorted to due to paucity of funds. 
The present study is aimed to help in precisely such situations 
when one has to judge between lining the sides only ox base 
only for partial lining of the channels, 

1 ,2 Why Numerical Method of Analysis of Seepage; 

Considering the importance of seepage through channels, 
considerable v/ork has been done by researchers to analyse the 
seepage through channels. Seepage depends on a vast number 
of factors, foremost being the soil type. It is difficult to 
accurately represent the actual soil strata occuring in nature 
in the theoretical analysis of seepage. Therefore, most of the 
earlier analytical solution and analyses , done by analogy 
studies and sand models, tended to simplify the actual physical 
problem so as to make it amenable to either mathematical 
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manipulations or easier experimentation. Such simplifications 
therefore provided solutions which v/ere merely estimates of 
the actual results. 

Numerical techniques, which are of recent origin, 
are comparatively more versatile and capable of handling 
complex physical conditions. Though actual physical situa- 
tion may be approximated to facilitated quicker results, they 
are essentially not simplified. Keeping in view this versati- 
lity of the technique and also the availability of the resour- 
ces for the numerical analysis, it was chosen for the present 
study , 

There are two powerful numerical analysis techniques 
available; namely, the Finite Difference Technique and the 
I inite Element (FE) Analysis. Both are quite versatile and 
are capable of providing quick and accurate solutions. However, 
one over riding advantage of the finite element analysis over 
the finite difference technique is its ability to quite accur- 
ately approximate the physical boundary of the region and the 
possibility of applying boundary condition to the whole 
system rather than to individual elements. Because of this 
advantage it is generally used for seepage analysis problems. 
Using the FE method, the solution of the basic equation 
governing two-dimensional steady state unconfined seepage 
through porous media is obtained for various sets of para- 
meters affecting seepage when either only the base is lined 
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or only the sides are lined. The results are discussed in 
chapter 5. 

1.3 Organisation of the Present Study ; 

The present v/ork is presented in six chapters and 
an appendix. Chapter 2 deals with the brief review of the 
past work on seepage through channels. The literature review 
is presented under three heads, namely, analytical studies, 
analogy studies and numerical studies aiongvjith the limita- 
tions of each . In the numerical study category reference 
has been made to most of the finite element analysis for 
unconfined flow through porous media. The chapter ends with 
an emphasis on the need for the present study. Chapter 3 
deals with the formulation of the governing equation for 
unconfined steady state seepage into the matrix form. The 
basic assumptions of the seepage through channels, Darcy’s 
law and the Laplace equation, brief description of the finite 
element method, the eventual formulation of different matrices 
and application of boundary conditions are described. Genera- 
tion of phreatic surface and calculation of actual seepage 
quantity are also discussed toward the end of the chapter. 
Chapter 4 is related to the description of the computer 
program ' CANSEP-86 ' used for the present study. Important 
subroutines of the program are discussed. An automatic Mesh 
Generation technique which is incorporated in the program is 
also described. This program was developed at I.I.T, Kanpur 
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by Achar [23], and certain necessary modifications were 
introduced for the use of the present study as described in 
chapter 4, Chapter 5 presents the results of the present 
study on seepage through partially lined channels. The 
verification of the program with the results of the previous 
workdone, the choice of the relaxation factor and optimum 
mesh size are discussed in the beginning of the chapter. The 
results of the present study are presented in tho form of 
graphs and tables and are discussed. Chapter 6 lists the 
conclusions from the present work and recommendations for 
further work on this topic. Appendix deals with the instru- 
ctions for giving the input for a given problem, two examples 
at the end further clarify the instructions. 
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CHAPTER 2 

LITERATURE SURVEY 


2.1 Introduction ; 


Steady state two dimensional flows through 
porous media which obey Darcy's Law are essentially governed 
by the Laplace equation. The problem of seepage through 
channels have been tackled by a number of investigators. Most 
of the solutions available are analytical. Analytical solutions 
for partially lined canals using transf ormation techniques or 
certain analogy methods are available too. However, most of 
the numberical studies available on seepage are related to 
earth dams. Nevertheless, the basic requirement in the pro- 
blems of seepage through earth dams or through channels laid 
in porous media is to obtain the phreatic surface which is 
not known a priori. Therefore, the formulations and approach 
to the problem in the formar case can be adopted for the 
present study, albeit with the required changes in the bound- 
ary conditions. A brief review of the literature available on 
various aspects of seepage analysis using different methods 
of analysis is presented in this chapter. 

2 , 2 Literature Survey : 

The previous work done on the subject can be 
broadly classified into three categories s- 



(i) Analytical studies, 

(ii) Analogy studies, 

(iii) Numerical studies. 
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The literature available within each category 
IS described below. 

2,2,1 Analytical Studies: 

Analytical solutions are available for evalua- 
ting seepage losses from channels under steady conditions for 
the following cases: 

i) Channels located in homogenous and isotropic 
medium extending upto infinite depth with shallow water table. 

ii) Channels located in homogenous and isotropic 
medium extending to finite depths with shallow water table. 

iii) Channels located in homogenous and isotropic 
medium extending upto infinite depth with very deep water table. 

Theoretical solutions for steady state seepage 
from unlined canals of various cross sections and for differ- 
ent types of flow boundaries are presented by various invest- 
igators, notably Vedernikov, Polubarinova - Kochina [ij, 

Harr [2], Muskat [3], and Hammad [4], All these solutions 
are based on highly simplified physical situtation which is 
quite amenable to mathematical manipulations. 

Garg and Chawla [6j have obtained a closed 
solution of seepage from a trapezoidal channel in homogenous, 
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isotropic media extending to infinite depth with drainage 
at a finite distance from the channel. 

Sharma and Chawla [7], using conformal mapping 
method gave a solution of seepage from channel in a homogenous 
medium extending to a finite depth to drains located at a 
finite distance from the channel considering horizontal and 
vertical drainage. In case of horizontal drainage, the stream- 
lines of seepage flow from channel join the bottom of the 
drain at various points along its width. In case of vertical 
drainage, seepage might enter from both sides of the drain. 

The streamline of seepage on both flanks of the drain would 
reach the drain in a horizontal direction. 

Lining of channels is now an accepted practice 
prevent 

to reduce seepage losses ^nd^water logging , Therefore partial 
or complete lining of channels deserve a detail study. Analy- 
tical solutions of some of the problems of partially lined 
trapezoidal channels using hodograph and Schwarff- Christoffel 
transformation techniques have been reported by Subramanya 
et al.[5]. The main cases studied are: i) Seepage from 
trapezoidal channels lined at sides only and bottom unlined, 
ii) Seepage from trapezoidal channels lined at bottom only 
and sides unlined. The medium is assumed to be homogenous, 
isotropic and of infinite depth. Also, the effect of water 
table is not taken into consideration and it is assumed that 
the capillary action is absent. The variation of the 
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seepage quantity with the aspect ratio and side slopes of 
the channels in the aforementioned two cases vias studied. 

It was shown that for shallow channels, the bottom lining 
v^as more effective than side lining and resulted in about 
50 percent reduction of seepage. For deeper channels it v'as 
seen that side lining v/as more effective. 

Limitations of Analytical Studies: 

The seepage quantity given by the theoretical 
formulae corresponds to steady state conditions which are 
difficult to attain in practice due to fluctuations of water 
depth in the channel and the water table. Thus ex cept for 
cases where steady state conditions can be attained early 
as for instance in soils of high permeability with main 
drainage located close to channel, the discharge as computed 
from the theoretical formulae may be poor estimate of the 
likely losses and only serve to indicate the order of seepage 
losses for deciding the necessity or otherwise of lining. 

The assumption of homogenous and isotropic soil 
in the complete seepage zone may not be always true. Als 
the boundary conditions created by presence of nearby pumping 
wells etc. render the application of analytical solutions 
trivial . 

2,2.2 Analogy Studios: 

For solutions of practical problems such as 
(i) seepage from a channel underlain by a highly pervious 
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matorial or impervious material at a finite depth from the 
base of the channel, (li) channel underlain by layered soil 
of different permeability and (iii) for the case of channel 
bed being practically impervious as a result of silting, the 
exact solutions are not available. In all such situations, 
an alternative to solv? the basic governing differential 
equations, namely, the Laplace equation with the usual 
boundary conditions, is by using analogy methods like electri~ 
cal analogy, sand models, electrolytic tank etc. 

Analogic studies of unlined and partially lined 
canals for both steady and transient flow conditions fusing 
resistance network analog have been reported by Bouwer [8,9j, 

A study was made of the theoretical effect of depth and shape 
of channel and position of the ground -water table on seepage 
from channels or streams. Soil conditions included in the 
analyses were uniform soil, soil underlain at varying depths 
by much more permeable soil, by much less permeable soil, 
or by a free draining layer, and a thin, slowly permeable 
(clogged) layer at the wetted perimeter of the channel. The 
contribution of unsaturated flow to seepage was examined 
and analysed. 

In actuality, seepage is a dynamic process that 
IS complicated by the non-uniformity of soil, by the effect 
of water quality, erosion, sedimentation, biological activity, 
etc,, on soil hydraulic conductivity, by fluctuating water 
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levels in the channel and in the ground, periodic dry ups 
of the channel , etc. For an analysis of how seepage is 
affected by conditions of channel, soil and groundwater 
table, etc. simplifications v/ere introduced and the following 
basic seepage models vjere distinguished ( Refer Fig. 2,1): 

Condition A~ Seepage into uniform soil underlain by 
material of very high hydraulic conductivity. If the water 
table is at or below the top of the permeable underlying 
material, the condition reduces to the case of seepage to 
drainage layer . The condition has received considerable 
attention in literature and is refered to as condition A'. 

Condition B ~ Seepage into uniform soil underlain by 
material of much lower hydraulic conductivity. 

It may be noted that the same seepage models have 
been considered for the present study. 

In the context of the present study it will be 
worthwhile to mention the important observations of Bouwerfll ,12] 
regarding the seepage through partially lined canals, though not 
all of them have been substantiated with actual experiments. 
According to Bouwer, lining of channel banks alone did not have 
a significant effect on the water table position, and thus 
not on seepage, in a viscous flow model study of seepage under 
condition B with relatively large [13 ] . For smaller values 
of D^, more significant seepage reductions can be expected. 
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Lining of tho channel bottom alone does not result in 
appreciable seepage reductions, except for small values of 
and/or Dp when most of the original seepage already takes 
place through the bottom. Also, the presence of free draining 
layers at relatively small depths below channel bottom can 
help tho effectiveness of tho bottom lining of a channel. 

However, there was no literature available which 
substantiated with experimental or theoretical data the 
important observations discussed above and therefore in the 
present study it has been attempted to establish the validity 
of fact. 

Other analogy studies in the related field are 
the study of seepage from leveed rivers into low lying 
adjoining lands by Todd and Bear [l4 ] with a tank electric 
analog and study of seepage through partially lined canals 
using electrolytic tank by L, R. I , P,R. I. [15 Amritsar, 

Limitations of Analogy Studies 

With analogs ,anisotropy is not readily handled 
and transformation to an equivalent isotropic medium is 
necessary. Thus, the complexity of the actual 
physical situation is compromised. 

Studies with sand models are particularly prone 
to scale errors and also equipment and experimental errors. 

Such studies are generally confined to two 


dimensional problems only. It is particularly difficult to 
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apply diverse boundrsry conditions when analysis of seepage 
is done using resistance analog technique. 

2.2.3 l^lume^lcal Studies; 

Flow through porous media have been analysed 
using the numerical methods like the finate element method 
and the finite difference method quite extensively. The 
works available relate to the problem of analysis of free 
surface flows . Since seepage through channels also comes 
under the same category the literature available on free 
surface flows is listed alongwith the literature related to 
seepage . 

Zienkiewicz and Cheung [l6] have presented a gen- 
eral method of analysis o flow through porous media using 
the finite element technique , Taylor and Brown [l7] have 
obtained solultion for flow through porous media with free 
surface. They used the fluid pressure as the dependent vari- 
able and the flow rate was determined after the pressures were 
obtained, An approximate free surface was assumed initially 
and the pressures at various nodes obtained. The pressure 
at nodes on the free surface was checked and if it was not 
equal to zero or greater than the permissible error, the 
whole procedure was repeated after modifying the position 
of free surface until the pressure on free surface became 
zero. Other works for location of free surface are presented 
by Liggett [l8] and Baiocchi, et al« [19]. 
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Finn [ 20 ] presents a finite element analysis of 
seepage through dams. The region of steady seepage was 
divided into contiguous triangular elements. The nodal poten- 
tial was used as the variable. The nodal potential on the 
free surface was compared v/ith the elevation of the node above 
a datum. The procedure was repeated until the potential of the 
free surface node became equal to its elevation. Askev/ and 
Thatcher [21 ] have presented an algorithm for calculating 
the seepage through a porous dam using the finite element 
method. 

Neuman and Witherspoon [22] performed an analysis 
of steady seepage with a free surface. Issacs [24] has 
reviewed various methods of adjusting the phreatic line in 
seepage analysis. It is concluded by him that the method used 
for the adjustment of phreatic line makes no significant diff- 
erence to the results obtained in principal application and the 
choice of method depends, therefore, on its suitability for the 
particular application and the ease with which it can be imple- 
mented , 


Das Gupta and Mustafa [31 ] have adapted the gen^-ral 
purpose program, FEAP, developed by Taylor [26] to convert 
the package into a continuous simulation for seepage problem 
with a free surface boundary, The results obtained are veri- 
fied with Bouwer’s [ll] analog study on seepage from channels. 
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Varoglu and Finn [27] present a variable domain 
finite element analysis of seepage from a triangular ditch 
into permeable soil underlain at a finite depth by a drain. 

A finite ^element analysis of seepage from a trape- 
zoidal channel towards a ditch has been reported by Bxebbia [28]. 
Two cases v/ere considered; a) permeable wetted perimeter and 
b) permeable sides with impermeable bed. The soil system was 
assumed to consist of a layer of isotropic, homogeneous top 
soil underlain by an anisotropic layer. The base under this 
layer was considered impermeable. Like Bouwer Brebbia 

too does not provide any data to substantiate the various 
observations regarding the reduction is seepage. 

A solution of transient free surface flow in porous 
media is given by Shun, Cheng and Li [29]. The numerical scheme 
is applied to isotropic and anisotropic earth dam problems and 
also to drainage from a ditch. 

Seepage flow across a discontinuity in hydraulic 
conductivity is presented by Issacs [30]« 

Limitation of Numerical Techniques 

Even the most efficient finite element computer codes 
require a relatively large amount of computer memory and time. 
Hence, use of the technique is limited to those who have 
access to relatively large, high-speed computers. 
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Tho results of a numorical rnc-thod must be interpreted 
with care and one must be a’vare of the assumptions employed 
in the formulation etc. 

The most tedious aspect of tho use of the finite 
element method is the basic process of subdividing the cont- 
inuum and of generating error - free imput data for computer. 
Although these processe^may be automated to a degree, they have 
not been totally accomplished by computer because some engin- 
eering judgement must be employed in the discretization. 

2 , 3 Need for Present Study : 

Loss of water through seepage from irrigation channels 
constitutes quite a substantial percentage of tho total utili- 
zable water. It has been estimated that transit losses in the 
alluvial channels in the r >rthen India are 17 percent for the 
main canals and branches, 8 percent for distributaries, and 20 
percent for the water courses which gives a total loss of 45 
percent of the water entering the canal head [38], The necess- 
ity of lining irrigation channels canot, therefore, be overemph- 
asised. Besides reducing seepage losses in channels and thereby 
leading to more efficient utilization, the li.-ing prevents 
waterlogging and weed growth and reduces requirement of land 
and maintenance cost. It has been estimated that lining of all 
the existing unlined channels in India could save enough water 
to irrigate an additional six million hectares against the 
present area of 36.4 million hectares [38], 
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Ho'vever, the complete lining of irrigation channel? 
involve huge costs. A practical idea would be to partly 
line the channel at least as a first phase, if comparative 
reduction is obtainable. As seepage from canals depends on 
a large number of variables, experimental or numerical analyses 
are generally preferred for obtaining seepage. As no data are 
available for seepage through partially lined channels, the 
present study is aimed to provide the same. Analysis has been 
performed using different values for the depth of v;ater table, 
depth of pervious or impervious layer, head in the channel and 
either base or side lined conditions. The data obtained are 
presented in graphical as v;ell as tabular form for easy ref- 
erence to people actually designing the lining of channel. 
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CHAPTER 3 

FORI/iULATIOH OF THE GOVERNING EQUATION 
3 . 1 Basic Assumptions : 

The basic assumptions that govern the flow through 
a porous medium are: 

1. The flow occurs only in the saturated region of a 
porous soil and that the pores in a medium such as sand 

r 

are interconnected and that the fluid can flow through 
them . 

2. The flow is laminar such that the Reynolds number of 
the flow, given by 

R = ■ ; where d = average pore diameter 

|i = fluid viscosity 
i' = fluid density 
V = velocity 

IS less than unity. 

3. Both the fluid and porousmedium are incompressible. 

4. The viscosity of fluid is constant. 

Other assumptions for the present study are: 
a) The flow is steady and two-dimensional i,e. the 
seepage zone extends upto infinity on either side 
of the section considered without any change in the 
material or boundary conditions. 
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b) The bottom and side of the channel are not 
clogged * 

3 , 2 Governing Equation; 


The relationship governing ground v^ater flow is 
known as Darcy's law. Its form can bo deduced from the 
general Navicr-Stokcs equations under the assumption that 
the flow is laminar and that the inertial forces are negligi- 
ble wlien compared with the viscous forces. Darcy's law 
establishes a linear relationship between the velocity of 
fluid flow and the hydraulic gradient. It may be stated as 



( 3 . 1 ) 


whore , 


0 = + h 


V = seepage velocity 

k = coefficient of permeability 
p = pressure 

Y = specific weight of water 
h = elevation above datum. 


In case of isotropic soil, a symmetric seepage tensor is 
introduced 


‘ V 1 
' ^ 1 

•^11 

ki 2 

•^13 


d0/dx 

- 

^21 

*^22 

^23 

1 


d0/dy 


' ^31 

^32 

k ^ 
*^33 J 

1 ' 

1 

1 1 

60/62 1 


If the principal axes of this tensor, xx, yy, zz are chosen 
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Eq. (3.2) bocom'-s 


V 

y 


V 


z 


0 0 ; d0/dx I 

; I 

0 k 0 b0/dy , 

yy \ . 

0 0 ' d0/dz ; 


(3.3) 


V/hen the seepage is concidcxod as tv^o-dimensional and the 
principal axes are horizontal (/-direction) and vertical 
(Y-direction) , then 



For homogenous media 


■XX 


>v 


( 3 . 7 ) 


If there is a source or sink of capacity Eq (3.7) can be 
written as 
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Equation (3.8) is the governing differential equation for the 
two dimensional steady state ground water flow through a 
hofiiogonous medium. 

By the application of the principles of variational 
calculus, it may be shown that the function, 0 , which mini- 
mises the functional, I, where, 



over the whole region of solution is the solution of the 
differential equation (3,7) . 

3 , 3 Finite Element Method ? 

For an exact solution the value of 0 would have 
to be determined at each and every point within the region. 

In the finite element method the aim is to determine 0 at a 
finite number of specified points or nodes. Between the node 
points, 0 is assumed to vary in such a way that it can be 
uniquely defined in terms of nodal valUes, Lines may be 
drawn between the node points which subdivide the whole 
region into a number of small regions or elements. Then the 
assumed variation of 0 over any one element can be specified 
in terms of the values at nodes within or along the boundary 
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of the element. In the program used here, the element 
chosen is triangular because any shape can be simply appro- 
ximated by triangles and graded meshes in which the distances 
between nodes vary depending on the accuracy and detail requ- 
ired are easily produced. Since there are only three nodes 
associated with each element , 0 is assumed to vary linearly 
within the element. The algebraic functions used to inter- 
polate the field variable over an element are called inter- 
polating functions. 

Often polynomials are selected as interpolating 
functions for the field variable in view of the ease with 
regard to their mathematical manipulations like differentia- 
tion and integration. The degree of the polynomial chosen 
depends on the number of nodes assigned to the element and 
the nature and number of unknowns at each node. The inter- 
polating function should be so chosen that the following 
general requirements are met: 

a) At the element interfaces or boundaries the 
field variable 0 and any of its partial derivaties upto one 
order less than the highest order derivative appearing in 
the functional must be continuous. This is called'compati- 
bility' requirement. 

b) All uniform states of 0 and its partial deri- 
vatives upto the highest order appearing in the functional 

(e) 

I (0) should have representation in 0 , when, in the limit, 
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the element size shrinks to zero. 


Considering any element (e) in the region 
D (Fig. 3,1 a) the function ^ v;ithin tl’^e element can be 
represented in terms of the nodal values as 

where is the variable within each element and 

represents a column matrix of the nodal values of the 

. r 0 ) 

function and J is the row vector of interpolation 

functions, sometimes known as shape functions also, which are 
functions of the coordinates of the nodes. The function 
defining the state variable v;ithin each element must satisfy 
the compatibility and completeness criteria as outlined above. 


The development of natural co-ordinates for tri- 
angular elements ( Fig. 3,1 b ) follow the choice of the 
co-ordinates also known as area co-ordinates, to 

describe the location of any point P(x ,y ) within the element 
or on its boundary as given by equation (3, 9) , 


The original cartesian co-ordinates of a point 
in the element should be linearly related to the new co-ordina- 
tes by the following equations. 


X = L.Xj^+LjXj+L|,X|^ 


(3.10) 


y = L. y. +L .y .+L. y. 

11 J J K K 


(3.11) 
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A third condition, rocjuiring thst tho v.'Gighting functions sum 
to unity, IS 

= 1 (3.12) 

so that only two of the natural co-ordinates can be independent. 

The point (^pfYp) within the triangle ijk divides 
the total area into throe seporate areas as shown in (Fig. 3, I b), 
The ratios of those areas to the total area define uniquely the 
position of P within the triangle. These area co-ordinates 
are : 


= A./A 

S = A/" 

‘-k = 


(3.13) 


These ar a co-ordinates ar related to th_ cartesian x and y 
co-ordinates by the following equations 


Li = (a^+bj^ X +c^ y)/ 2 A j 

L. = (Sj+b^ X +c^ y)/ 2 A (3.14) 


where. 

^i 





^i 

= Yj-Yk 

bj^Yk-Yi 

(3.15) 


c . 

1 


Cj=X.-Xk 



(3.15) 
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and 


or 


2A 


equation 



(3.16) 



’'k Vk 

can be ro-v/ritten as 


(3.i6a) 


= [( ^ ^ (3,16b) 

The field variable 0 may be interpreted as a function of L , 

Lj , instead of x,y and z as 

0 ^®^ = 

It may be noted that equals 1 at node i and equals 0 at 
nodes j and k. and varies linearly between those points. 


Also , 


60 ^ 
6x 


(g) 


■5a ( h'ii + * \«>k) 


ilL 

dy 


(g) 


(3,18) 


= sr <"A + "k?*k> 

: function 0^®^ is thus uniquely and continuously defined 
throughout the region in the terms of nodal 0 values and there - 
fore also satisfies the compatibility and co mpleteness requ- 
irements , 


3.4 Formulation of the Governing Equation in Matrix form; 

The governing differential equation for steady two 
dimensional seepage flow, as mentioned earlier is 
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H) ^ 37 If) + 3 =0 (3.19) 

with the boundary conditions 

0 = 0(x, y ) on Sj^ , and 

^ ''"^y fy '^y '^2 the usual notations. 

For the location of phreatic surface, the conditions on the 
phreatic surface that are to be satisfied are 


and 


i) 0 


li ) 


dn 


h { i.e. p=0) v/here h = elevation head 

0 i.e. the flow across the phreatic surface 
is zero. 


of (i) and (ii ) only one of them is used as a boundary condi- 
tion to obtain the solution of phreatic surface from the FE 
equations and the other is used as a check of the solution. 
Depending upon this check the iterations are performed on the 
equation. 


If 0 is represented as equation (3.9) over an 
element, then the residual from equation (3. 19) , is 


R 


(e) = (k M 


(e) 


•) + 


X - “57 

According to Galerkin's criterion, 


(K 


^1^) + Q (3.20) 


fl 

^(e) 


dx dy = 0 


(3.21) 




FIG, 3.2 EXTREME TOUNDARIES 


y. 


dy=- B Cos 0 
=-clB 



dx=d8 Sin© 

=d8 Oy 


FIG. 3,3 DIR3C'SI(n'‘AL COMFONfci.'TS 
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Substituting for R' in the above, we get 

■ ■''i [57 C'x '> * + QJ dx dy = 0 


dy "y <^Y 


(3.22) 


Let Xj^ , X^, Y^, and Yj. represent the extreme boundaries of 
the domain D, in the x and y directions respectively as shown 
in Fig. 3,2. 

On integration by parts, Eg, (3.22) becomes 


rr N. [ ^ (k 
*^(6) ^ ^ 


“d— + 5y dr“) + Q J dx dy 


= ■/' (“1 ''x ’ '“y - -f-f, ''x 


dx dy 


+ J (N. k 


Y. 

^ I * 


y 57 ‘''i''' ''y "V" 


+ // 0 dx dy 


(3.23) 


From Fig. 3.3, we have dy = ~ dB n, 

and dx = dB n 


Applying Stokes’ transformation, 


(e) , ^r 


and, 


f N. 

\ " ax Xl 


;(e) 


dy = / N.,B kjj dB 


(3.24) 


{ = /VBV 

^ b 


(3.25) 
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Substituting equations (3,24) and (3,25) in Eq. (3.23) and 
regrouping, we get 


= If Q n dx dy + / tl , B (k, 


— n 

X X 




ff (k 'N,x> + k |N,v' ) dx dy 

^\e) ' X ’ ■ «' V ' ' ay 


=// Q .N, dx dy + / ’N.Bi (k n 

D(e) 


(3.26) 


For the triangular olement, 0'®'^ is presented as 
0^^^ = a + bx+cy 


then 




(3.27) 

(3.28) 




(3.29) 


Equations (3.27) to (3.29) show that there is completeness 
of n^*^ “1 order and compatibility of order, 

1 

//(k^ N,X|- [N,xJ + k .N,yj [N,yJ ) dx dy |0 } 

D(e) 

=// n:N|dx dy -/0;N,B|qdB (3.30 


(3.30) 



Condition A { Condtt 
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or 


C I = , Q - / q ;j,D d3 


(3.31) 


whoro , 


[k 0 ^®^] = SI (k^ r,x; [U,x] + k H,y; [ri,y]) dx dy 

) 


D 


' k b, b +k c c . 
X 1 1 y 1 1 


'i° ^ k b b 4-k c c 


! k b b. +k c c, 

; X 1 k y i k 


where , 


k b b +k c c 
X i j y 1 j 

k b b -f-k c c , 
X 3 3 y 3 3 

k, b,b, +k c c, 
X ] k y J k 


k b b| -t-k c c, i 
X 1 k y 1 k , 

i 

k b, b -f-k c, c . I 
X k ] y k 3 

I 

k b, bi +k c, c, I 
X k k y k k 


(3.32) 


bi=yj-y^ ’ 






=k= d"’*! 


3 . 5 Boundary Conditions : 


For a complete specification of any problem, 
certain boundary conditions must be given. There are gener- 
ally two cases: 

(a) the value of 0 is specified, This presents no 
difficulties and the value of 0 is prescribed 
and kept fixed at the appropriate nodes. In the 
present study, the boundary of the region of 
analysis with the GWT, the boundary of unlined 
part of the channel and the bottom layer if fully 
permeable constitute this boundary condition 
( Fig. 3.4). 
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The negative term, /q ' M, B > dB , on the right 
hand side of Eq. (3,31) n^nifests this boundary condition. 
If q = 0, this term vanishes and the boundary condition that 
the velocity vector is tangent to the external boundary is 
automatically satisfied by the formulation above, 


3 , 6 Generation of Phreatic Su rface : 

Before the ea"ation (3.31 ) is finally solved, 
using the Gauss elimination technique, the boundary conditions 
discussed above are applied to the global stiffness matrix that 
had been assembled for the complete domain. The solution of 
Eq, (3,31) gives the 0 at each node in the domain. Though 
ZienkiGV'acz [32] and others have reported that, within the 
numerical accuracy, the finite element method gives the correct 
value of 0 throughout the region for fixed boundaries, in the 
present case the upper boundary i.e. the phreatic surface is 
not fixed and therefore the boundary conditions on it may not 
be at once fulfilled. Therefore, either a trial and error 
method, reported by Finn [ 20j or an iterative method, as 




PI6.3.5 C0NVER6ENCE CRITERION 
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adopted in the present study, may be made use of. 

In the iterative method an initial trial position 
IS first chosen by the user and the FE analysis done. The 
results from that analysis can be used to generate automati- 
cally the next position of phreatic surface and the cycle 
repeated until a solution within numerical accuracy is ach- 
ieved. It can be seen that this method is far less tedious 
on the part of the user. However, to automatically generate 
the phreatic surface nodes, the complete mesh too has to be 
regenerated. To this end, an automatic mesh generation tech- 
nique was incorporated by Achar [23] in the program which was 
used for the present study after some modifications. The 
algorithm for automatic mesh regeneration is the same as 
given by Durocher and Gasper [25], The test of convergence 
for various relaxation factors (a) was carried out by Achar[23j . 
The result is reproduced in Fig 3.5 and has been adopted with- 
out change. Relayjation factor {a) is the factor applied to 
the differenCv,- between the 0 of the phreatic surface node 
obtained after FE analysis and the y-coordinate of that node 
on the mesh. The now position of the node is obtained by 
changing the old value by the product obtained above 

or + a .6 ; whoro, (3_34) 

a— Relaxation factor 
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Calculation of Velocity and Di _sc harqe : 

Though emohaois is more on obtaining 0 (i.e, 
piezometric head ) at -/arious points in a domain using the 
IE technique, for the present study the velocity of flow and 
the actual quantity of discharge (or seepage) across a 
boundary are more important points of interest. The computa- 
tion of these quantities from the primary variable at the 
node i.e. 0 is as follov/s: 


By Darcy's la'w, 

v= - k 1 


(3,35) 


writing Eq. (3.35) in the matrix notations, for a two dimen- 
sional flow, 


, v 


X , 


:k 0 ! ;d0/dx i 


I 


t : Vyl “ ' 0 ky' ',d0/dy > 

But the gradient matrix is written as 

bk : : 01 ; 

Ck : ; 0, 


(3.36) 


I 60/ dx , ' b, 

/ ' 1 ' ^ 

! d0/dy!“ 2A I 






(3.37) 


01 


where i>jfk represent the three nodes of the triangular 


element and 0^ etc. the heads at node i etc. A is the area 


of the element 
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Substituting c'^u?tion (3.37) in equa tion (3. 36 ) , we get 


V 

X 

1 


0 : 


‘’k ; i i^i 


> 

l( 

> 

2A 

-0 

k . 

y 

' c . c 

1 J 

i 

(3.38) 








The quantity 

of flow 

IS 

given 

by 




< 

> 

11 

a 







= k 1 

A 




(3.39) 

The 

quanti ty 

of 

f low 

can also 

be related to 

the 


seepage velocities bv computing equivalent flows at the 
nodes. For this purpose, the net flovjs in the x and y dire- 


ctions 

at the nodes 

are 

evaluated as 

the uniform 

flow tributary 

to the 

node and we 

get 





' Q • 
j ; 


f 

0 



) xj / 


i ^ 

1 

0 i 



Qxk 1 


! 

0 ; 


Iq! = 

: ; 

i 

= 

1 1 0 

2 

i 'J V ' 

1 J 

(3.40) 




i 0 

1 

1 \ 

bi 

1 



1 Qyk I 


! 0 
> 

^k 



where b.,b ,b, and c. ,c.,c. are constants for a particular 

1 J K 1 j K 

element and .Vf is as calculated in Eq. (3.37). 

If 


























FIG. 4.2 FLOV, CHART FOR B-SE LIHED CASE 










FIG. 4.3 FLOIV CHART FOR A* CONDITION 
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CHAPTER 4 
THE PROGRAf-1 

The computer program ' CANSEP*-86 ' used for the 
present study was developed by Achar [23], It is a versatile 
program for dealing v/ith the two dimensional steady state 
groundwater flov/, Non-homogeniety and anisotropy can be 
easily handled. For the purpose of present study,arnodif i- 
cation was made to incorporate the condition of seepage 
through channels lined at base only in which case there is 
a possibility of a bottom phreatic surface being generated. 
Also, to deal with A' condition (Fig. 2.1) of seepage where 
the phreatic surface meets orthogonally the underlying 
permeable layer certain changes v/ere introduced. The phreatic 
surface in such a condition requires not only vertical adjus- 
tment so that he head of the node on the phreatic line is 
equal to the y co-ordinate of the node, but also horizontal 
displacement so that the condition of orthogonality is ful- 
filled, The flow chart for the main program is given in 
Fig, 4.1. The flow charts of the changes made are given in 
Fig. 4.2 and Fig. 4,3. The actual method of giving data 
input is explained in Appendix I with two examples. 

4.1 Important Subroutines; 

The program consists of twelve subroutines out 

of which eight subroutines are actually used for the present 
study. 
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Following is the brief description of the various subroutines 
used . 

i) Subroutine DATA: 

This IS the most important subroutine as it 
reads the complete specification of the problem, any error 
in giving the input can either lead to spurious results or 
abortion of the execution of the program. The geometry of 
the region of analysis, the various boundary conditions, 
material properties, the total number of nodes and elements, 
number of free surface nodes, etc. are all required to be 
specified. This is the only subroutine in the complete 
program that reads the input data, therefore any changes to 
be made in the style or content of input may be made in this 
subroutine only without searching the whole program for READ 
statements . 

li) Subroutine MESH; 

The user specifies only the co-ordinates of the 
nodes of the superelements into whach the region of analysis 
is initially divided. This subroutine subdivides the super- 
elements into specified number of nodes, number them and 
codifies their x and y co-ordinates. The actual formulation 
of elements from these nodes is done by subroutine LINEAR* 
Since this subroutine is derived from the work presented by 
Durocher and Gasper [25], a brief description of Automatic 
Mesh Generation Technique is given below: 
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Zienkievvicz and Phillips have described the 

theory and technique employed in a mesh generator based on 
isoparametric mapping concepts 4 A typical parabolic super- 
element is shov^n in Fig. 4, 4, the x and y co-ordinates of a 
typical point within the superelement is related to the 
eight pairs of nodal co-ordinates by the equations: 


n 


E N. (p.p) 

1=1 ^ ^ 

(4.1) 

n 

E N (p,Tj) y 

1=1 ^ ^ 

(4.2) 


where, (|3,h) is the shape function associated with node 
'i’. The shape functions for parabolic quadrilaterals as well 
as other element types have been given by Zienkiewicz [32], 
Thus, if the nodal co-ordinates (x^,y^) are known, the carte- 
sian co-ordinates of any specified point (p,ri) can be easily 
calculated from Eqs. (4.1) and (4,2). In case of the three 
noded triangle used for the present study, the eight noded 
superelement shown in Fig. 4,4 suffices. The relatively simple 
programming logic can be used to generate a mesh of a given 
refinement once the user specifies the number of element 
subdivisions in the p and t) directions. 

ill) Subroutine LINEAR: 

This subroutine divides the noded region of a 
superelement into finite elements which are triangles in the 
present study. It takes into consideration the fact that 
equilateral triangles produce more accurate results than 
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long narrow triangles by comparing the two diagonals of the 
quadrilateral produced by four nodes by subroutine WESH. 

It also generates the element connectivity array and material 
specification codes. 

iv) Subroutine HEAD: 

It generates the element stiffness matrix using 
the equations obtained in chapter 3. The element stiffness matr 
ices are then assembled into the global stiffness matrix. 

By the help of subroutine SIMUL, the boundary conditions are 
applied to these equations. The subroutine SOLVE then solves 
by Gaussian elimination these simultaneous set of equations to 
obtain the nodal heads. 

v) Subroutine SliViUL: 

This subroutine introduces the boundary condition. 
The head wherever specified as the boundary condition is 
multiplied by the corresponding stiffness to obtain the 
nodal discharge. This discharge is then added to the right 
hand side of the equation reproduced below. 

|k} L0J<''®> = LQJ 

Thus the total number of simultaneous equations is reduced to 
the number of nodes with specified Q only and thus the set of 
equations can be solved to obtain 0 at different remaining 


nodes . 
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VI ) Subroutino SOLVE: 

This subroutine solves the* matrix equation 
obtained by subroutine SIMUL using the Gaussian olinnination 
method . 

vii) Subroutine VELO: 

This subroutine calculates the gradient (d0/ds) 
over the element from v;hich the velocity of flow over the 
element is obtained using Eqs. (3.37) and (3.38). From the 
elemental velocity the nodal discharges are calculated using 
equation (3.40). To obtain the total seepage from the channel 
the nodal discharges of the nodes lying on the boundary of the 
seepage zone are added up. For the present case, since the 
symmetry of the seepage zone is exploited by considering only 
half part of the zone fox analysis, the total seepage quantity 
obtained by this subroutino is also half the actual seepage, 
Therefore, care should be taken when the total seepage quantity 
is required. 

viii) Subroutine OUTPUT: 

This subroutine is called after the iterations 
have been completed and the phreatic surface(s) located. All 
the necessary information like element connectivity, heads 
at each node, nodal discharges, elemental velocities, co-ordi- 
natesof the phreatic surface, total seepage, etc, is printed out 
by it. Like the subroutine DATA, this also is self -contained 
and any change in the output may be made in it without 
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searching the whole program for '.VKITE statements. 

4 . 2 Some Physical Specifications of the Program: 

i) The program can handle 200 nodes and 200 
elements at a time. This limit is only due to the dimension 
allotted to the array of the nodes and elements. Larger 
numbers can be accommodated by increasing the dimensions of 
the arrays. 

ii) The rnesh generation program limited to the 
main program can generate a wide variety of elements ranging 
from 3-noded linear triangles to 8-noded isoparameteric 
quadrilaterals. However, it is important to note that the 
program can be used for analysis by using 3-noded linear 
triangles only because the subroutine HEAD incorporates 
the stiffness matrix bast-i on the interpolating function 
for linear triangles only. Therefore, necessary changes in 
subroutine HEAD will have to be made to accommodate higher 
and complex elements. 

iii) The total CPU time needed to obtain the 
phreatic surface depends largely on the initial position 
estimated by the user. Thus, a proper engineering judgem- 
ent must be exercised at the time of specifying the phreatic 
surface initially so that it is as close as possible to the 
actual position. However, it has been observed during the 
present study that the convergence becomes increasingly 
difficult with deeper GVVT and also when GVtfT is at or belov/ 
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the underlying permeable layer and CPU time as high as tv;o 
minutes with approximately ninety iterations was required 
at times in such cases. The average CPU time was 30 seconds. 
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CHAPTER 5 

SEEPAGE FROM PARTIALLY LIHED CHANNELS 

5 . 1 Introduc tion ; 

Using the FE program ’ CAMSEP-86' described earlier, 
the characteristics of seepage from unlined and partially lined 
channels under different boundary conditions were studied. The 
details of this investigation and the salient characteristics 
of seepage from partially lined channels are described in this 
chapter. 

5.2 Basic Flow Conditions; 

In actuality, seepage flow systems are characterized 
by channels of irregular cross section, non-uniformity of soils 
in horizontal as well as /ertical extent, changing elevations 
of the water surface in the channel and the water table in the 
soil, and other complications. Perfectly impermeable lining 
of channel sides or bottom is an exception rather than the 
rule. Therefore, in view of the multitude of complexities 
with which the seepage flows are afflicted in nature, theoreti- 
cal treatment must begin with the approximation of the soil 
and boundary conditions. The accuracy of the seepage and 
seepage reduction predicted for a given channel in that way 
then depends largely on how well the pertinent soil, water 
table, and boundary conditions can be characterized. Although 
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the FE analysis used for the present study yields accurate 
answers, the fact that the model for which the solution is 
obtained is always a simplified version of the field situation, 
renders the solution an estimate at best. 

5.2.1 Theoretical Models: 

There are two basic conditions to v/hich the multitude 
of natural profiles of soil hydraulic conductivity can be re- 
duced for theoretical treatment of seepage flovj systems. These 
conditions labeled A and B in accordance with Bouwer [ll,12], 
are as follows : { Refer Fig. 5,1) 

Condition A: The soil in which the channel is laid is 
uniform and underlain by more permeable (considered here as 
infinitely permeable) material. 

Condition B : The soil in v/hich the channel is imbedded 

is uniform and underlain by less permeable (considered as 
impermeable) material. 

The case of seepage to a free-draining permeable 
layer in the subsoil is a special case of condition A, and it 
is obtained by letting the water table be at or below the top 
of the permeable material. This condition is labeled A*. 

5.2.2 Significant Parameters Governing Seepage: 

To obtain general solution, the depth of infinitely 
permeable (Dp) or impermeable (D^^) material for condition 
A and B respectively is treated as a variable. The position 
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TABLE 5.1 

Range of Parameters Used in the Present Study 


Variable Description Par'^moter Range 


Hw 

Head of Water in the 


0.25,0.50,1.0,2.0 

channel 


Dw 

Position of vyatcr 

table from the water 

surface in the channel 


1 .0 to 3.0 
depending 
upon HyWj^ 

Dl 

Position of the under- 
lying impermeable 
layer from the base 
of the channel 


1 .0,2.0, 3.0 


Position of the 
underlying infinitely 
permeable layer from 

the base of the 

channel 

V'-’b 

1 .0,2.0, 3.0 




m 

Side slope of 

the channel 

m 

1.0 
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of the v^ater table at some distance from the channel in also 

treated as a variable (D ). This water table acts as a solid 

v<i 

boundary for the present case of steady state seepage. Other 
variable is the head of water in the channel. All these 

parameters are non-dimcnsionali zed with refcienco to the width 
of the base (Wj^) of tho channel. Thus, for a given channel 
( say a trapezoidal channel of side slope, m ) the seepage 
quantity Q is written as, 

Q = fn (H^,D^,D^ or Dp, , k ) . Refer Fig. 5.1 for 

the definitions of these variables. Non-dimensionalization yields 
n I. K, D D, 

VJ Vk ~ ^ ^ WT' VT ’ W~ ^ ^ ^ (5.1) 

s b b b D 

The value of m is kept constant for the present study as unity. 
The value of k (=0.01 m/s) is also kept constant. The non- 
dimensionalized parameters and their ranges considered for the 
present study are given in Table 5.1. 

5.3 Description of Problem; 

Thus having simplified the actual conditions of seepa- 
ge into two basic theoretical models and defined the governing 
parameters and their ranges of variation, the problem for this 
study reduces to obtaining the absolute seepage quantity for 
each of the models for all the possible combinations of the 
different values that each parameter can take for the channel 
being unlined, base only lined and sides only lined cases. 

Having obtained the absolute seepages, the percentage 
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reduction in seepage due to base or sides lining of the 
channel is calculated as, 


Percentage reduction in seepage = 
where , 


— — X 100 / 


0 . 


UL 


(5.2) 


^UL ~ seepage from unlined channel 

Qj^ = seepage from lined ( base or sides) channel. 


For a given case of lining and a set of parametric 
conditions, the actual seepage is obtained in two stages. In 
the first stage, using the formulations describe in chapter 3, 
the phreatic surface is determined. And in the second stage, 
the seepage quantity is determined by summing up the nodal 
discharges across a given section, 

5.4 Determination of the Phreatic Surfac e; 


In determining the phreatic surface, the heads of 
the nodes supposed to lie on the phreatic surface are compared 
with their corresponding y co-ordinates. This is to satisfy 
the boundary condition of zero pressure head on the phreatic 
surface of an unconfined ground water flow system. If the 
difference between the y co-ordinate and the head is greater 
than 0.01 meter then the y co-ordinate is modified using the 
expression 

>'new = ''old + (^'node-V'old* • “ <^-3) 



Ynev; ~ niodafied y co-ordinate. 

Void ~ previous y co-ordinate, 

^node~ node on tlie phreatic surface, 

a = relaxation factor. 

The cycle of nev/ mesh generation obtaining the head of all the 
nodes in the region of analysis and the comparison of heads 
and y co-ordinates of phreatic surface nodes and the modifica- 
tion of y co-ordinates of nodes lying on the phreatic surface 
using Eq. (‘'.3) is repeated until for all the nodes on the 
nodes on the phreatic surface the difference betv/een the 
y co-ordinates and the heads of nodes on phreatic surface is 
less than 0,01 meter. 

For the case when the bottom of the channel is also 
lined, there is a possibility of the generation of the bottom 
phreatic surface. In such condition, the solution is said to hav 
been obtained only when the condition of difference between the 
heads of the nodes on the phreatic surface and their y co-ordi- 
nates being less than 0.01 meter is satisfied simultaneously 
for both the general phreatic surface and the bottom phreatic 
surface. 

5.4,1 Chojce of the Relaxation Factor; 

The choice of the relaxation factor governs the 
convergence of the problem. A poor choice may cause large 
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oscillations of the values obtained about the actual value. 

For the present study, the o obtained by Achar [23] after 
conducting a test { Fig, 3,5) was adopted without modifica- 
tion, after having been checked and found satisfactory 
keeping in viov; the accuracy achieved and the CPU time recjuired. 
The a adopted for the present study is 0.25, The samea is 
seen to hold good for the bottom phreatic surface also, 

5.5 Seepage Quantity; 

The equations for obtaining the quantity of seepage 
in matrix form are given in the section 3,7 of chapter 3, The 
total quantity of seepage across the channel is obtained by 
summing up the nodal discharges in both x and y direction 
across any selected section. As for the present study the 
symmetry of channel geometry and surrounding soil and water 
table conditions is exploited thereby permitting analysis of 
only half of the region of seepage, the discharge so calculated 
is also half of the actual quantity of seepage taking place 
through the channel. Therefore, if absolute quantity of 
seepage through the channel is requiredi then the quantity 
obtained from the program ’CANSEP -86' should be multiplied 
by two , 

5.6 Optimum Mesh Size: 

The size of the mesh plays a vital role in the accur- 
acy of the results obtained by the FE analysis. It is an 



Percentage of error 


^0 



Number of elements 


F10 5.2 TEST FOR OPTIMUM NUMBER OF ELEMENTS 






FI6 s.3 verificat:on of results of present study for 

UKLINED CHANNEL WITH BOUWER^ (lt,12 J RESULTS 
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present study is not to stipulate precisely the amount of 
seepage reduction by lining the sides or the bottom of the 
channel but to provide a close estimate of the actual redu- 
ction in seepage that may take place. 

5 . 7 Verification of the Progr am; 

It v/as essential to establish the authenticity of 
the program before it could be put to further use by ^/erify- 
ing the results produced by it v/ith some established v.'ork 
done earlier, Bouwer [ll,12] has presented a set of curves 
for seepage through unlined channels using electric analog 
method for various values of the parameter Bouwer had| 

in turn, verified his results with the results obtained by 
investigators like Dachler, Ernst, Vedernikov and Hammad [4], 
Therefore, Bouwer's results can be assumed to be reliable and 
verification with them can be taken as the proof of the authen- 
ticity of the program used for the present study. 

Fig. 5,3 shows a comparison of the values obtained by 
using the program and Bouwer's [ll,12] values. These curves 
are for an unlined channel of side slope 1H:1V and =0.5, 

The variable is plotted on the x-axis while the depen- 

dent variable Ig/k (Eg. 5,1) is on the y-axis. The curve 
parameter is condition B and A respectively. 

The values are in good agreement with Bouwer's results, However, 
it may be noted that when the values of the variable are 

such that the water table is quite close to the underlying 



SeepoQtf in percent of origin a 



H^/Wb 


FIG. S4. VERIFICATION OF RESULTS OF PRESENT STUDY FOR 

CHANNEl LINED AT THE BOTTOM WITH BOUWER^ 
(11 ,12 ) RESULTS 
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bottom layer (Refer Fig» 5.1) the difference betv;een 
Bouwer*s results and the results obtained by the present 
study increases. The reason for this may be that v/hen the 
position of the groundwater table (GVVT) belov' the channel is 
quite deep the hydraulic gradient becomes steeper but as the 
mesh size remains constant for all positions of GV/T the 
percentage error increases slightly above normal. The diff- 
erence is of the order of about 3/, 

Fig. 5.4 shows a comparison of the reduction in 
seepage obtained by lining the bottom of the channel v/hen the 
underlying bottom layer is highly permeable (condition A) for 
the parameters indicated as obtained by Bouwer [ll,12] and 
by using the program. Again, the results are in excellent 
agreement, 

5.8 Seepage Through Partially Lined Channels; 

The program ' CANSEP-86’ was executed repeatedly for 
unlined, sides lined and base lined cases for a given set 
of parameters indicated in Table 5,1 to obtaine the seepage 
quantity in each case from which the percentage reduction in 
seepage was obtained using Eg. 5.2. The conditions A and B 
were treated separately. 

5.8.1 Seepage for Channels with Sides Lined and Base Unlined; 

Seepages for the channel with sides lined for the 
complete range of and for condition B and of and 
for condition A were obtained. The side lining of the channel 
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FI6. 5.5 TYPICAL MESH SHOWING SUPERELEMENTS AND 

FINITE ELEMENTS 





FI6 5.6 TYPICAL FLOW SYSTEMS OBTAINED FOR SEEPAGE 

IHROUGH CHANNEL WITH SIDES LINED AND BOTTOH 
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was achieved by specifying zero nodal discharges at nodes 
3,22 and 23 of the FE mesh ( Fig. d, 5). The condition B was 
simulated by specifying zero nodal discharges at nodes lying 
on the bottom layer while condition a necessitated specifying 
the heads equal to the GWT at the nodes lying on tfie bottom 
layer. The position where the phreatic surface met the GV/T 
was taken as ten times the width of the base of the channel[l2j, 
Since the results are plotted in terms of percentage reduction 
in seepage with respect to the seepage through uniined channel 
the seepage quantity for the same set of parameters and condi~ 
tions as for the side lined case was obtained for the unlined 
case. Typical flow systemsfor condition B, condition A and 
condition A' are shovm in Fig. 5.6. The results of the analysis 
are presented in terms of dimensionless graph? showing percentage 
reduction as a function of for different values of 

(condition B) and Dp/Wj^ (condition A) (Fig. 5. 7 a,b,c). The 
curve joining the end points of the curves obtained for condi- 
tion A and different values of Dp/Wj^ represents the condi’ion A' 
when the GWT lies at or below the underlying infinitely perme- 
able layer. It may be noted for A' condition, any 

further lowering of water table does not increase the effective 

value of D . Thus, for condition A', the D , values at the 
w w 

abscissa should be interpreted as D^+H . The analysis was 

p W 

performed for trapezoidal channels with side slope of 1H:1V 
and four different water depths (expressed as The 

values are presented in Table 5.2. 
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The curves, in Fig. 5.7, for porcontago rf.duction in 
seepage for a channel with sides only lined shoy; that for 
condition B there is very little effect of D on the percen- 
tage reduction of seepage. However, for concJirion A, 

V/ 

significantly affects percentage roauction in seepage but 

than 6 times Wj^ only. Therefore, for condition A, 
the effect of water table deeper than 6 times can be ignored 
without causing appreciable error in the percentage reduction 
in seepage obtained due to side lining of the channel. 

The curves indicate the depth of impermeable layer 
from the base of the channel (D^) for condition B tends to 
increase the percentage reduction in seepage v/ith increase in 
the . However, this increase in percentage reduction in seepa- 
ge is very less. For condition A, the depth of permeable layer 
from the base of the chan' il (D^) plays a significant role in 
the percentage reduction in seepage. The increase in causes 
significant decrease in the percentage reduction in seepage. 
However, the curves for percentage reduction in seepage versus 

the D /Wl- tend to come closer to each other as the D„/\Vu 
w b P u 

values increase. As can be seen in Fig. 5.7 a, the curve for 

D /K =3 IS almost a straight lino, indicating thereby that 
p b 

any further increase in would not significantly affect the 
percentage reduction in seepage, Similar trends are apparent 
for curves in Fig. 5,7 b and Fig. 5.7 c» 

The percentage reduction in seepage for a given set 
of conditions increases with the values when the sides of 
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the channel are lined as can be seen from Fig. 5.7 a,b and c, 
5.8,2 Seepage for Channels with Base Lined and Sides UnlinerJ; 

Seepages for the channel with base lined for the 
complete range of D and D for condition B and of D and D 

1 v/ P v/ 

for condition A were obtained. The base lining of the channel 
was achieved by specifying zero nodal discharges at nodes 1,2 
and 3 of the FE mesh (Fig. 5,5) . The conditions A and B were 
simulated as for the channels with sides lined, The position 
where the phreatic surface met the GV/T was taken as ton times 
the width of the base of the channel. The results are presented 
in the form of non-dimensional parameters same as those used 
for the sides lined channels. Typical flow systems for condi- 
tions B, A and A' are shov/n in Fig. 5,8. The percentage 
reduction in seepage obtained for each of the set of parameters 
and condition are given in Table 5.3. The results are oresented 
in the form of graphs in Fig. 5.9 a,b,c. Condition A* represents 

D =D + H , Further lowering of D., does not increase the effe- 
w p w VJ 

ctive value of D , Thus, for condition A' the D values at the 

W Vv 

abscissa should be interpreted as Dp+H^, 

The percentage reduction in seepage due to lining of 
the base for condition B very small as seen from Fig. 5, 9 
and Table 5,3 as compared to the percentage reductions for 
condition A. The effect of the position of water table 
is marginal on percentage reduction in seepage for condition B. 
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Thus, the effect of water tabic pocition on th. percentage 
reduction in seepage can be neglected fox condition B. 
However, the effect of on percentage reduction of seepage 
for condition A is significant and the percentage reduction 
tends to increase Vv'ith increase in D /W. . 

W D 

The depth of impermeable layer (D^) for condition B 
affects the percentage reduction in seepage by increasing 
it when itself is also increased. On the other hand, 

an increase in Dp/Wj^ for condition A tends to decrease the 
percentage reduction in seepage. 

5 . 9 Instances of the Development of Phreatic Surface at 
the Base of the Channel ; 

When the base of the channel is lined with an imper- 
meable material, there is a possibility of the flow line 
seperating from the base and developing into a phreatic 
surface. This phenomenon, however, is expected to depend 
upon the permeability of th- medium in which the channel is 
laid(k), the width of the base the position of the GV/T 

(D ) and the position of impermeable (D ) or permeable layer 

(D ) underlying the channel- During the present study it was 
P 

observed that for condition B, there was no instance v;hen the 
phreatic surface was developed. This fact was further explored 
by increasing Wj^ and and decreasing k but still the base 
phreatic surface did not develop* On the other hand, for 
condition A, for the normal working dimensions and values ^ 



of parameters > same as those used for condition B, the 
base phreatic surface promptly developed. 
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Therefore, conclusion can be drawn that though one 
might expect the development of beise phreatic surface v^hen 
the base of a channel is lined, its generation is significan- 
tly affected by the permeability of the underlying layer of 
soil . 

5.10 Alternate V/ay of Representation of Results ; 

Though partial lining of the channels leads to 
some amount of reduction in seepage, the reduction is not 
proportional to the percentage of the wetted perimeter of 
channel lined. That is, if, say, 60 percent of the channel 
perimeter is lined it is not necessary that reduction in 
seepage should also be 60 percent. Often it is less. This 
fact IS highlighted by obtaining a factor called the lining 
factor (fL)* The lining factor is the ratio of the change 
in seepage quantity per unit unlined vjatted perimeter of cha- 
nnel to the seepage quantity per unit wotted perimeter of 
unlined channel. Thus, 

Q^/Pt - Ql/(Pt" (5.4) 

Qy/PT 

= Lining factor 

_ Seepage quantity for unlinod channel 


where, 




Hw/Wb 

I6.S.10 percentage reduction in seepage by lining the 

SIDES AND BOTTOM SEPARATELY FOR AN UNDERLYING 
PERMEABLE LAYER 
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~ Seepage quaneity for partially lined channel 
“ Total wetted pori''nL;tcr of unlincd channel 
~ Total wetted unlincd portion of channel. 

For a linear relationship butvyu-n th>. perccritage reduction in 
seepage and the percentage of the lined portion of the channel, 
the value of f^ should be zero, hov.ever, if the percentage 
reduction is less than the proportion of lining than the 
value of f^ would be negative. Thus, grater than zero 
represents a cost efficient channel while f^ less than zero 
represents a cost inefficient channel, Table 5 . 2and 5.3 represent 
the lining factors for all set of parameters and conditions 
and cases of lining obtained by this study. 

5,11 Comparative Performance of the Tv;o Types of Partial 
Li ninqs : 

From the results of the present study presented in 
section 5.8 for the base lined and sides lined cases, it is 
apparent that there is substantial reduction in seepage 
(ranging from 30-55 percent) due to lining of either base 
or sides only for condition A as compared to the reductions 
in seepage for condition B { ranging from 2-15 percent) for 
the range of parameters considered. Owing to this fact, the 
curves represented in Fig. 5.10 are prepared for the partial 
lining cases for condition A only. 

The curves of Fig. 5.10 represent the percentage 
reduction in seepage due to partial lining of the channel 
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versus the depth of water in the channel (represented non- 
dimensionally as for various values of Por 

each set of curve parameters, the curves of sides lined only 
and base lined only cases are obtained. These curves inter- 
sect at one point for each set of parameters. The point of 
intersection is very crucial for deciding the depth of v/ator 
in the channel such that for greater than the value for 

this depth, sides lining is more effective than base lining 
and for less than this value base lining is more effe- 

ctive than sides lining in reducing seepage. From the curves 
of Fig. 5.10, a value of about 0.40 for H /W. is obtained 

Q w' b 

for Dp/W,^ varying from 2 to 3 and pp = 2 such that for 
less than 0.40 the lining of the base only would bo more 
effective than lining the sides. The percentage reduction 
achievable is greater than 33/. depending upon the H /VL value. 
However for greater than 0,40 sides lining would be more 

effective in reducing seepage and the percentage reduction is 
greater than 33;^ depending upon the value of H /Vi . 

W D 

Since the curves shown in Fig. 5.10 demarcate the 
limits of efficiency of the base lined and sides lined cases 
of partial lining of a channel, thoir importance cannot be 
overemphasized while deciding the type of partial lining for 
a channel. For a given set of conditions of Dp, and V/j^, 

the curves help to obtain the most economical solution subject 
to the condition of a given percentage reduction in seepage. 
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CHAPTER 6 

CONCLUSION AND RECOM^/iENDATIONS 


6,1 Conclusion ; 


With transit losses m irrigation channels 
upto 45 percent of the total intake at the canal head, the 
importance of lining of the channels cannot be overemphasized. 
Though considerable research has been done on analysis of 
seepage through unlined channels, there was no literature 
available on analysis of seepage through partially lined 
channels except certain observations reported by Bouw'er [11,12] 
and Brebbia [23], Thus a need was felt for further study in 
this area. 

The analysis of seepage through partially lined 
channels for a set of different variables that govern seepage 
was done using a program ' CANSEP-86', initially developed 
by Achar [23] but was modified for present use, based on finite 
element analysis. Before being put to use for the present 
study the results obtained by the program were verified with 
works of Bouwer [ll,12] on seepage through unlined channels. 

Steady state two-dimensional unconfined seepage 
through channels depends upon a host of factors and though 
all such variables and various boundary conditions can be 
taken care of by the finite element method of analysis, the 
actual physical situation was approximated by two conditions 
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A and B keeping in vievv tho general nature of thii study. 

The side slope of IH : IV v/as kept constant in the present 
study. The following significant conclusions based on the 
results of analysis performed for base lined only and sides 
lined only cases for each of conditions A and B arc obtained: 

1. It is observed that the condition A is 
more active than condition B and therefore larger quantity 
of seepage as well as higher reduction in seepage is obser- 
ved for condition A for the cases of both base lined only and 
sides lined only. 

2. There is very marginal effect of the position 
of water table (D ) on reduction in seepage through base 

W r 

for 

lined only as well as side only lined channels^condi tion B and 
therefore can as well be reglected while analysing seepage 
reduction. 

3. The percentage reduction in seepage decreased 

with increase in tho depth of infinitely permeable layer 

(D ) for condition A while increased with increase in the 
P 

depth of impermeable layer (D^) for condition B, and, as 
expected, the results show a trend of approaching a common 
value of percentage reduction in seepage for both conditions 
A and B when and approach infinity. 

4. The percentage reduction in seepage is not 
found to be directly proportional to the percentage of the 
wetted perimeter lined for partially lined channels. For 
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deeper channels (H^/\‘}^=2) lining factor as lo^; as - 2 has 
been observed for sides only lined cases for condition A 
thereby severely affecting their cost efficiency though 
sides lining only may reduce seepage to the order of 60 
percent . 

5. The development of bottom phreatic surface 
was not observed for any of the cases for condition 6 when 
the base only was lined while it was readily apparent for 
condition A, thereby indicating the role played by the 
permeability of the underlying layer. 

6. For condition A» for Dp/V/j^ =2 and 
and for Dp/Wj^=3 and Dyw^j=2 it is observed that for 

less than 0.4 base lining only of the channel is more effective 

in reducing the seepage than side lining only. However, for 

H greater than 0,4 sides lining only is more effective in 
w b 

reducing seepage than base lining only > The order of magni~ 
tude of percentage reduction in seepage for 

may be as high as 70-80 percent for base lining only case when 
H /W,_/lGss than O. 06 and upto 55 percent for the case of 
sides only lined for condition A when is greater than 2. 

6.2 Recommendations ; 

1, As has been shown by the present study, the 
set of curves obtained for different cases of lining, namely, 
base only lined and sides only lined .representing the 
relationship between the percentage reduction in seepage and 
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very helpful in deciding the limits of the effici- 
encies of each of the lining cases for a set of different 
values of parameters in terms of It would be v/orthwhile 

to obtain more such sets of curves for large ranges of values 
of the parameters governing seepage so as to be useful for 
practical purposes, 

2, The side slopes of the channel wore kept 
constant for the present study» The analysis may also 

be performed for various values of side slopes which are 
generally provided like 1,5H:1V, 2hl:lV etc. 
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APPENDIX I 

PROGRA// USER IHSTRUCTIOMS 
The following assumptions were made for this 

program. 

1. The positive direction of x is from loft to 
right and y is postive upv/ards with the underlying layer , 
permeable or impermeable, as the datum. 

2. The principal axes of seepage coincide with the 
X and y axes. 

3. The units for length and time are in meters and 
seconds respectively. The pressure is in terms of the 
head in meters of water above the bottom layer, while the 
discharge is in cumecs. 

All the input data are format free. Numbers 
must be separated by one or more blank spaces or by a comma. 
The program is written in FORTRAN -IV. 

The general instructions are described below. Two 
different examples given at the end further clarify the 
procedure, 

I, General User Instructions; 

Set of 

Input Data 

Data 

1 (i) INRG: Total number of superelements in which 

the area of analysis is divided. 

(ii) INBP: Total number of superelement nodes. 
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1 (ctd, ) (vii) IREMUfA: Renumbering of thj node numbers and 

the elements. 

= 1 Renumber the mesh. 

= 0 Do not renumber the mesh.. 


(viii) IBOTLM: '/hether or not the chrnnel is 
lined at the bottom? 

= 1 Bottom of the channel is not lined. 

= other ottom is lined, 

(ix) IDIFIM: Vfliether there are different 
impermeabilities in the area of analysis? 

= 1 Soil is homogenous. 

= other soil is not homogenous, 

(x) lADASH: l/hether water table is at or belov/ 
the permeable layer? 

= 1 V/ater table is not at or below the 
permeable layer. 

= other V/ater table is at or below the 
permeable layer. 


2 ( To be read if IPOL is not equal to zero) 

(i) IPOLAR: Polar co-ordinate reference number, 

(ii) XOFSET (IPOLAR): Offset of the origin of the 
polar co-ordinate system in x-direction with 
reference to the global co-ordinate's origin. 

(ill) YOFSET (IPOLAR): offset in the y-direction 

of the origin of the polar co-ordinate system 
with reference to the global co-ordinate's 
origin. 
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(i) IPOLAR; Co-ordin^ite reference. 

— 0 Cartesian co— orainato s/stem. 

— 1 Polar co-ordinate system reference 

number 1. 

■ 

i 

= M Polar co-ordinate system reference 
number N. 

(ii) K: Superelement node number. {K=1,INBP) 

(iii) XP{Kj1)j X CO— ordinate of the node if 
IPOLAR = 0 else radius of the node position 
from the origin at the datum, 

(iv) YP(K,1): Y co-ordinate of the node if IP0LAR=0 
else angle of the node position from the 
origin at the datum. 

4 ( To be read only if ITMP is not equal to zero) 

(i) K: Superelement node number (FC=1,INBP) 

(ii) XP(K,2): The nodal temperature at the super- 
element nodes. 

5 ( To be read only if ITHK. is not equal to zero) 

(i) K: Superelement node number. (K=1,INBP) 

(ii) YP (K,2): The nodal thickness, 

6 ( To generate the elements) 

(i) NRG; Number of the superelement, 

(ii) NRWI(NRG); Number of divisions in the 
superelement columnwise. 
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6 (ctd. ) 


7 

8 


(iii) NCOLI(NhG); Number of aivisions in the 
superelement rovmiso, 

(iv) INDl : Indicates the total number of boundary 
nodes for the given superelement, 

= 1 Total numbv^r of boundary nodes 
equal to eight, 

= 2 Total number of boundary nodes 
equal to four. 

(v) N0P(NRG,K),K=1,8,IND1 : The superelement 
node numbers read consecutively and in 
counter-clockv;ise direction. 

(vi) NMAT(NRG): Material code number for this region, 

(i) TITLE: Any title for specifying the problem 

(maximum 16 characters), 

(i) UNIT: Units in which various quantities 
are represented. 

=. 1 M.K.S. System. 

= 2 F.P.S. System, 

(i) TYPE: Type of the media, 

= I Homogenous and isotropic porous media. 

= 2 Non-homogenous and anisotropic porous 
media. 

= 3 I'llhen the phreatic surface traverses 
porous media of different hydraulic 

conductivities , 
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10 

(i) 

NNl; Total number of nodes of thr finite 
element mesh. 


(ii) 

NE : Total number of elements into which 

the region is divided. 


(iii) 

NST: Total number of nodes at which thio 

head is specified. 


(iv) 

NSQ: Total number of nodes at which flux is 

specified. 

11 

(i) 

QQQ: Flow generated in the element initially. 

12 

(i) 

NT: Number of the node v/hlch has head 

specified. 


(ii) 

TNT: The head which is specified at a certain 


node (nx). 

(NT and TNT are read in pairs till the total 
number of pairs is equal to MST'', 

13 ( To be read only if IMPDIF is not 'equal to 1 ) 

(i) TKU: Permeability in the upper (right) 
region of the area of analysis. 

(ii) TKD.' Permeability in the lower (left) 
region of the area of analysis, 

14 (i) LIMITL: The node number such that all node 

numbers lesser than it have permeability 
TKU and greater than it have permeability TKD, 


(ii) LIMITU: The node number such that all node 
numbers greater than it have permeability 
TKU and lesser than it have permeability TKD, 
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15 (To bo bead if IMPDIF is equal to 1} 

(i) TKX (l), 1= 1,NE: Permeability in x-direction • 

(li) TKY (I), 1= lyNE: Permeability in y-direction 

NFS: Total number of phreatic surface nodes. 
NFY (I), 1= 1,NFS: Node numbers, related to 
the initial super elements, of the nodes 
lying on the phreatic surface, 

(ii) NFT (l), 1= i,NFS: Node numbers, corresponding 
to the finite element mesh, of the nodes 
lying on the phreatic surface. 

18 (To be read if IBOTLN is not equal to 1 ) 


19 


(i) NOBNOD: Total number of nodes lying on the 
bottom of the channel. 

(ii) DWATBL: Depth of the GWT from the datum, 

(i) NBPHY (I),I= 1,^^DB^J0D: Node numbers, 

r 

related to initial superelements, of the 
nodes lying on the bottom of the channel. 

(ii) NBPHT (I), 1= 1,NOB{'JOD: Node numbers, 
corresponding to the finite element mesh, 
of the nodes lying on the bottom of the 


channel. 
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